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Abstract. We discuss a numerical method to study electron transport in mesoscopic devices 
out of equilibrium. The method is based on the solution of operator equations of motion, 
using efficient Chebyshev time propagation techniques. Its peculiar feature is the propagation 
of operators backwards in time. In this way the resource consumption scales linearly with the 
number of states used to represent the system. This allows us to calculate the current for 
non-interacting electrons in large one-, two- and three-dimensional lead-device configurations 
with time-dependent voltages or potentials. We discuss the technical aspects of the method and 
present results for an electron pump device and a disordered system, where we find transient 
behaviour that exists for a very long time and may be accessible to experiments. 



Electron transport through mesoscopic devices contacted to leads is intensely studied in 
chemistry and physics (see e.g. Refs. [1,2]). The conceptual basis for theoretical studies is the 
non-equilibrium Green function (Keldysh) formalism [3]. The Meir-Wingreen- formula [4] allows 
for the calculation of steady state currents. For time-dependent potentials or gate voltages 
the calculation of a current is a serious problem already for non- interacting electrons. One 
crucial point in most approaches is that the resource consumption scales quadratically with 
the number of system sites, since each Green function Gij{t,t') or operator product c\{t)cj{t') 
has two site indices. The necessity to study large systems conflicts with the rapid growth of 
computational demands, especially if long leads with long recurrence time are required. Recent 
studies therefore addressed mainly one-dimensional (ID) situations, e.g. by time-propagation of 
operator expectation values or of several single-particle eigenstates [5-7]. 

In the present contribution we numerically solve the equation of motion for a single fermion 
operator. In contrast to previous studies, the initial time coordinate is propagated backwards 
in time. Then the computational effort, in particular the memory consumption, scales linearly 
with the number of lattice sites used to represent device and leads. This allows us to calculate 
the non-equilibrium current even for large systems that could not be treated otherwise. 

We consider the following situation: A device of x Ly x L^-sites is contacted to two long 
leads extending along the x-direction (see figured]). This slab geometry includes the 2D (L^ = 1) 
and ID (Ly = = 1) case. For the kinetic energy in the Hamiltonian 

= ^mCn + ^ J^4cm + Yl ^m(t)cLc„, " ^ ^Lc^ (1) 

^ ' x<l l<x<Lx x>Lx 

we assume nearest-neighbour hopping (oc t) along each axis i = x,y, z, and set t = 1 as the 
unit of energy. The potential energy terms describe the voltage bias U (t) applied between the 




Figure 1. Typical 2D slab geometry, 
witii a central device region contacted to 
a left /right lead. 

left /right lead, and the local potentials Vm{t) in the device. The current along a bond in the 
x-direction is given by the expectation value of the corresponding operator 

To specify the initial conditions, we assume that for times t < the device is isolated and 
uncharged, and the leads are in equilibrium at zero temperature (this situation can be described 
by setting i = along the contacts at x = l,Lx and letting Vm — > oo in H{t < 0)). At time 
t = device and leads are brought into contact, and electrons can flow onto the device. We 
then ask for the current I(t) = X],y2(V'(*)li(2; ^ 2) through a section x = const, of the slab. 

Time propagation of the many-fermion Fock state \ip{t)) is not feasible, and the reduction of 
this problem to the propagation of several single-fermion states is restricted to non-interacting 
electrons and introduces additional errors that must be controlled. It appears to be more natural 
to propagate the current operators themselves. With the time evolution operator U {t, t'), defined 
hyi§-^U{t,t') = H{t)U{t,t'), U{t,t) = 1, operators can be transformed to the Heisenberg picture 

A"{t,to) = U{to,t)AU{t,to) , (3) 
and expectation values are obtained from A^{t,0) and the initial state \il^o) = \^JJ{t = 0)) as 

mmm) = mA^'it^oMo) ■ (4) 

Therefore we must represent {t, 0) in terms of operators Aj = Aj^{0, 0) at time t = 0, 

A''{t,0) = ^aij,t,0)Aj , (5) 

j 

which allows evaluation of equation Q since all initial expectation values (?/'o|^j|V'o) are known. 
We observe that two ways exist to obtain A^{t,0). The standard way starts from 

i^^^(i,0) = U{0,t) [A,H{t)]U{t,0) = [A,H{t)f{t,0) . (6) 

In this equation, the commutator is given as a combination of operators Aj, which is evaluated at 
time t in the Heisenberg picture. Therefore solution of equation ([6]) requires time propagation 
of all operators Aj^{t,0). The total number of coefficients in the corresponding expansions 
equation ([5]) for these operators grows quadratically with the number of sites. 

To avoid quadratic growth we here proceed the opposite way, starting from —i-^lJ{t,t') = 
U{t,t')H{t'). In contrast to the previous case, we do not need to transform the commutator in 

^-^A^M = [H{t'),Uit',t)AU{t,t')] = [H{t'),A^{t,t')] (7) 
to the Heisenberg picture. Instead, it is only the single operator A whose expansion 



A''{t,t') = Y,<j,t,t')Aj 



(8) 



Figure 2. (Colour online) Current I{t) 
measured 10 sites to the left (black curve) 
or right (red curve) of a ID pumping 
device with Lx = 500 sites. The system 
parameters are V = 2, lo = tt/W, k = 
7r/2, with Fermi momentum kp = vr/4 
for the leads. Time t is given in multiples 
of the period T = lixjuj = 20. Shown 
is the net current I{t) = j^_rpl{t')dt' 
over one period, in comparison to the 
situation with a static potential wave 
(dashed curves). The insets display I{t) 
over the first and last 5 periods. 

must be propagated in time. We do not need to keep track of all the other operators Aj[t,t'). 
Therefore memory consumption and computational effort are proportional to the number of 
sites. Making use of (AB)^ {t,t') = {t,t')B^ {t,t'), each fermion operator in equation ([2]) can 
be propagated separately to obtain (jm)^(*i*')- 

The price to be paid here is that equation ([7]) leads to propagation backwards in time, 
evolving the second time coordinate from the initial A{t,t) = ^4(0,0) = ^ to the final A{t,0). 
The main reason why we nevertheless consider this procedure is that, avoiding quadratic growth 
of memory consumption, only this way allows for the treatment of large systems, which are 
inaccessible in the first way. With modern Chebyshev techniques [8], time propagation itself 
is fast and accurate, and the additional effort for backward propagation is compensated by 
the reduction of the problem size. Note that for periodic time dependence, operators can be 
propagated iteratively by a full period. 

Turning to the slab geometry, the eigenstates of decoupled infinite leads are discrete along 
the y-, z-direction and continuous along the x-direction. In the actual calculation, infinite leads 
are replaced by long leads of length L ^ Lx,y,z- The recurrence time revealing the artificial 
discretization is of the order L, and must be larger than the maximal time t in the calculations. 
Using operators Cq for lead eigenstates to energy Eq, and operators Cm for device sites m as the 
Aj in equation ^ the initial conditions specify the expectation values at time t = 0, 

(cLcn) = (cLca) = , (^C^/) = Saa'^ifJ- - e) , (9) 

where /i is the Fermi energy. Note that only these conditions change for finite temperatures. 
In our first example in figure EJ a travelling potential wave 

Vxit) = V cos{kx - oot) (10) 

pumps electrons in a ID system with zero voltage bias U = from the left to right (i.e. along 
the positive x-direction). In the beginning electrons flow onto the initially uncharged device. 
The initial transient behaviour evolves into a 'pseudo' steady state (left panel), which persists 
over 50 periods. The true steady state, with equal current at both sites of the device, is 
approached only in the extreme long time limit (right panel) . Comparison to the case of a static 
potential, with the same k but a; = 0, shows that the 'pseudo' steady state persists on much 
longer time scales than the initial transient behaviour. This may allow for the experimental 
observation of deviations from the steady state, which here occur over 300 periods up to the 
maximal propagation time. To allow for such long propagation time in our calculation, we 
choose long leads with L > 300 x T = 6000. The full system has + 2L > 12500 sites. 
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Figure 3. (Colour online) Current I{t) 
at the contact surface of a disordered 3D 
device with Ly = Lz = 3 and different L^, 
for 7 = 2 and kp = t^/^- A static voltage 
bias Ul = U/2, Ur = -U/2 is applied, 
with [7 = 1 (upper panel) and C/ = 0.1 
(lower panel). The inset compares the 
La; = 5 to the Lx = 50 system, and shows 
the time averaged current (dashed curve) 
in addition. Note that the data are given 
for a single disorder configuration. 



In our second example in figure [3l a constant voltage bias U transports electrons through 
a disordered 3D device, with (static) random potentials with uniform probability distribution 
£ [~7)7]- A steady state is approached only for the short system = 5 with large bias 
[7=1. Longer systems {Lx > 10) support localized states, whose presence causes oscillations 
in I{t). These become more pronounced for smaller bias U = 0.1. The net current through 
the device is close to zero, since localized states do not contribute to transport. Note that the 
average current for L^; = 5 is much smaller for [7 = 0.1 than for [7 = 1, while for the = 50 
system it is of the same magnitude. 

In the given examples, up to 18450 lattice sites are treated for device and leads. Because 
of linear scaling in our method we must keep only that many elements in memory, and all 
calculations can be performed on standard desktop computers. With quadratic scaling, we would 
have to store about ~ 1.7 x 10^ elements, requiring already 2.5 GByte for complex numbers in 
double precision. 

In conclusion, the numerical method discussed here allows for the calculation of time- 
dependent currents in large systems of non-interacting electrons, which are inaccessible to 
most other existing techniques. The propagation of operators backwards in time is a tolerable 
disadvantage in such cases. Apart from these achievements, the development of advanced 
solution techniques for complicated non-equilibrium Green function equations of motion appears 
more promising for future progress, especially with respect to the inclusion of dissipation or 
interaction. How such developments are possible within the context of the Chebyshev Space 
method [9] and Sparse Polynomial Space approach [10] will be discussed elsewhere. Data 
obtained with the present method serve as a reference to validate these new techniques. 
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